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BASIC IDENTIFICATIONS 

• Differential Equation System 

• Turbulence Closure 

• Numerical SolutIon Algorithm 

• Verification Tests 


ANALYSIS SPECIFICATIONS 


• Solution Domain 

• Boundary Conditions 

• Initial Conditions 

• Closure Model Constants 
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DIFFERENTIAL EQUATIONS 

CONTINUITY/ MOMENTUM AND ENERGY 

• First Order Effects 

Downstream Diffusion Negligible 
Transverse Momenta Governs Pressure 
Continuity Governs Transverse Velocities 

» Reynolds Stress Components of Equal Order 

• Well - Posed Problem Statement Required 

TURBULENCE CLOSURE 

• TKE - Epsilon Two-Equation System 

• Reynolds Stress Constitutive Model 
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TURBULENT THREE-DIMENSIONAL AERODYNAMIC FLOW 


3D Thin Layer Navier^Stokes Equations (1 $ L ^ 3) (2 £ & ^ 3) 
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three-dimensional directed flowfields - 

PARABOLIC ALGORITHM STRUCTURES 


Continuity ; 

o Pseudo - Pressure Correction 

Spalding et. al. (1970 - 78) 

9 Potential Function / Axial Vorticity 
Briley et. al. (1979 - 80), 

Ghia (1979 - 76), Dodge (1976 - 80) 

0 Penalty Differential Constraint 
Baker et. al. (1979 - 81) 

Pressure: 

• Axial Mass - Conserving Gradient 

"Partially - Parabolic" 
o Inviscid Poisson Equation 
Wall Modifications 

• Ordered Poisson Equation 

Complementary Solution (3D Potential Flow) 
Particular Solution (Turbulence + Convection) 
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TWO - DIMENSIONAL PARABOLIC NAVIER - STOKES 
STEADY, UNIDIRECTIONAL TURBULENT FLOWS 

FINITE ELEMENT PENALTY ALGORITHM 
Formulational Structure From Classical Variational Concepts : 
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Weighted Residuals Galerkin Algorithm 


Tip =Uii e (M) = U 


{N k b)f{(Ul 


J e 



L s (v h ] ’AU{M K }lifeJ') = {o] 


R 


R 


^Nk] L 1 « , e h , v h ) = {0} 


SS8290S00T 


Source: https://www.industrydocuments.ucsf.edu/docs/mjmkOOOO 



TWO - DIMENSIONAL PARABOLIC NAVIER - STOKES 
STEADY, UNIDIRECTIONAL TURBULENT FLOWS 

FINITE ELEMENT PENALTY ALGORITHM 


Penalty Function Form : 


LMeo) = V*u h * 0 

. 


Boundary Conditions On (|)^ : 

• Porous - Dirichlet 
« Non - Porous - Neumann 
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3DPNS FINITE ELEMENT PENALTY CONSTRAINT ALGORITII 
Finite Element Solution Algorithm: 
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3D PARABOLIC NAVIER-STOKES EQUATIONS 


SPACE MARCHING : 

• Uj, k, e, Ti 

• U 2 , U 3 ( Constrained S.Ti V » IT = 0 ) 

POISSON EQUATIONS : • 


• Pressure ( P c + P p ) 

• Continuity Function 


ALGEBRAIC 


• Equation of State 

• Reynolds Stresses 
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Yields 16 Variables t Node Point 
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VERIFICATION TESTS 


CONFIGURATIONS : 

• Bounded Flows 

Turbulence Closure 
Vortex Structures 

e Free Jets 

Entrainment 
Discreteness 
Vortex Structures 

• Shear Layers • 

Decay Rate 
Reynolds Stresses 



SOLUTIONS : 

• C' M C: 3DPNS Computer Program 

• CDC Cyber / 203 


6S8290S00T 

i 

i 

1 ; • 

Source: https://www.industrydocuments.ucsf.edu/docs/mjmkOOOO 




3-DIMENSIONAL PARABOLIC NAVIER-STOKES 

TURBULENT FLOW IN RECTANGULAR DUCT 

TRANSVERSE PLANE VELOCITY DISTRIBUTIONS 
(Baker, et. al., ASHE, 1981) 


HASS CONSERVATION ONLY 


EDDY VISCOSITY MODEL, Xl/D = 30. 



U m = 0.0003 
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THREE-DIMENSIONAL PARABOLIC NAVIER STOKES 
TURBULENT FLOW IN RECTANGULAR DUCT 

PROBLEM DEFINITIONS 



I98890S00T 
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3-DIMENSIONAL PARABOLIC NAVIER-STOKES 
TURBULENT FLOW IN RECTANGULAR DUCT 
TRANSVERSE PLANE VELOCITY DISTRIBUTIONS 


EXPERIMENTAL DATA, X/D ■ 37 

MELLING & WHITELAW (JFM, 1976) 
U M = .008 



FINER DISCRETIZATION M ° 25 z 
BAKER et a 1. (A l A A 1983) U 11 =.00A 
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Figure It, Penalty Algorithm Prediction of Turbulent 
Kinetic Energy Distributions, Four Multiple Jet 
Geometry, uj = 12 m/s, a) xj/R * 0.0, k° = 0.005, b) xj/U 
= 0.12, k m = 0.02, c)x|/R = 0,25, k m = 0.035, d) x|/R = 
1.5, k m = 0.039. 
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Figure 12. Summary of Penalty Algorithm Prediction 
of Axial Mean Velocity U| Decay as Function of Jet 
Initial Turbulent Kinetic Energy Level. 


Figure 14. Penalty Algorithm Prediction of Transverse 
Plane Mean Velocity U| Distributions as Function of 
Locator Radius Extension, Four Multiple Jet Geometry^ 
xi\ s 12 m/s, a) Gu^Milliine±c/ Extension^cj7 R - 0. 2 5, 
u'P = 0.1 S2, b) Extension/^*j7R = 0.7J, 

un* = 0.052, c) Extension, xj/R = 1.5, 

* 0.055, v — _ 
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Figure 15, Penalty Algorithm Prediction of Transverse 
Plane Velocity Distribution On Half Discretization of 
Full Solution Domain, uf = 12 m/a, x|/R = 1.5, uTC = 
0.084, M s 19 x 19. K 
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figure 16. Penalty Algorithm Prediction of Transverse 
*Hithe Velocity uj hiatr Ibiitioh, Pour Multiple jet 
Geometry with One 3et Off, Of = 12 m/s, xi/li = 1.3, 
tl'g = 0.126, 
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Figure 17. Penalty Algorithm Prediction of Smoke 
Visualization Distributions On Full Solution Domain, 

12 m/s, x|/R = 6.0, a) Four Jets Operating, Y m = 53%, b) 
Three Jets Operating, Y m * 39%. 
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THREE-DIMENSIONAL PARABOLIC NAVjER-SluKES 

MULTIPLE F R E E-J E T SECONDARY VORTEX FI 

3DPNS ANALYSIS DESCRIPTION 


Geometry / Discretization 


Transverse Plane Vortex Field 



' »»ii • i i ■•*••* 1 


'• —•! I I I * • i I I I 


* % i l ik. 


...all It# 


•Ml III I I I III* 

\ 1 11*1.. 9| j , , , , , , 9 * * 

S\W,. . ___ _ 



.i l I I i i l III*. 

.m ill i i I in. 




•**M I * < 


.•.il III • • • I I I® •«*•*» it i 


, .till III I III I I « | I 


O&8Z90S0OT 

I 

' I 


ELD 


Source: https://www.industrydocuments.ucsf.edu/docs/mjmkOOOO 












T&8290S00T 


I 

t 

I I 

Source: https://www.industrydocuments.ucsf.edu/docs/mjmkOOOO 



A e T R 0 N FILTER 


GENERATION 


OF VORTEX 'VELOCITY' FIELD 
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actron filter 

generation of vortex velocity field 
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THREE-DIMENSIONAL PARABOLIC NAVIER-STOKES 

multiple free-jet s e e o nd a r y vortex field 

EXTREMUM VORTEX VELOCITY FIELD COMPARISONS 


Smoke Flow Visualization 3DPNS Prediction 

i _ _ —- - --- 
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THREE-DIMENSIONAL PARABOLIC NAVIER-SWES 

, 1ULT ifle free-jet secondary vortex field 

initial turbulence uncertainty - X / R - !•? 

k 0 * 0.0025, u J 


0.005 , u £ 


influence of 

= 0.067 


k 0 » 0. U « = 0.009, 
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ACTRON FILTER 


EFFECT OF TURBULENCE ON AERODYNAMICS 


Vent Velocity Decay 


Vortex Velocity Maximum 
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actron filter 

VORTEX VELOCITY FIELD, X/R = 150X 
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ACTRON FILTER 

AERODYNAMICS MODIFICATIONS BY SLEEVED VENT 



Downstream X / R ( X ) 


_ Essential Action 

Nearfield 

• Enhance Vent Decay 

9 Prevent Entrainment 

q Strengthen Vortex 

Farther Field 

9 Dissipate Vent 

q No Entrainment 

o Weaken Vortex 
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AERODYNAMICS MODIFICATIONS BY SLEEVED VENT 


Vortex Velocity Maximum . Sensory Response Experiment 
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PREDICTION OF SECONDARY VORTEX FLOWFIELDS 
INDUCED BY; MULTIPLE FREE-3ETS 
ISSUING IN CLOSE PROXIMIT Y 

A. 3. Baker •, 3. A. Ori«chovski** t and G. E. Stungis* 

let ! ! 


A continuity Constraint Unite element numerical 
algorithm, for solution ol the three-dimensional 
parabolic NaviefrStokes equations lor subsonic turbulent 
Hows, is applied to prediction ol secondary vortex 
.llowfields induced by multiple tree-^w issuing in close 
•proximity. The combined action of axial decay ol the 
;high speed jets* induced entrainment from the farfield, 
♦ and geometric discreteness ol the jets yields prediction 
| pi eight counter-rotating vortex pair* lor a four-jet 
!Configuration. The extremum magnitude of the induced 
.transverse vortex velocity can exceed 10% ol that ol the 
j^t initial axial velocity. The results ol the numerical 
!predictions compare well with available video-graphic 
data, including the effects ol geometric and initial 
: flowiieid modifications. 


; t. Introduction ; 

• An important aerodynamics problem class is 
assessment ol the steady-flow interaction ol multiple 
subsonic free-jfts issued in close proximity. Dependent 

•upon geometric parameters, and whether the jets are 
laminar or turbulent* this uniaxial initial configuration 
Can induce detailed three-dimensional velocity fields 
characterized by large scale vortex structures. A 

• substantial challenge in computational fluid mechanics is 
to identify, and evaluate a numerical solution algorithm 
for prediction of this type of three-dimensional turbulent 
flow field. 

Figure 1 illustrates the essential geometry of a 
’representative aerodynamic multiset system. Provided 
there is no reversal of the dominant component of 
velocity* parallel to the axis of the jets, the three- 
dimensional parabolic Navies Stokes (PNS) equation 
system is a candidate for solution of the problem class. 
A formal order of magnitude PNS analysis, confirms! 
that axial diffusion processes are negligible, and that the 
.transverse momentum equations govern transverse plane 
pressure distributions. The continuity equation governs 
first order effects on momentum. For the Uee-jet 
'problem class, the complementary solution to the 
.resultant pressure Poisson equation is a homogeneous 
;constant, where the Poisson equation is obtained from 
I the divergence of the transverse momentum equations. 

• A turbulence closure model is alio required, and the 
^parabolic form of the two-equation \ turbulence kinetic 
.energy-isotropic dissipation system is representative of 
i the minimum acceptable level of simplicity. 


boundary conditions. A primitive variables formulation 

* must rearrange the continuity equation to yield a 
' deterministic system for transverse velocities, cf. Biker 

et.al.!, Dodged, Patankar A - Alternatively, a vector 
potential function can be defined to identically satisfy the 
. continuity equation/ and a vorticity equation derived to 
replace the transverse momentum equations, cl: Bhiley 
et.al. 4 , Ghia et.allV One dominant factor controlling this 
Basic decision is that the multiple free-jet solution 
domain is unbounded. Hence, llowfield specification in 
. the farfield is a priori unknown, ie-, the entrainment 
j action of! the jet interaction that induces inflow/efflux 
: across surfaces located arbitrary distances from the jet is 
a solution outpuu Fori this required flexibility, and the 
I inherent mathematical robustness, the finite element 
I penalty, function algorithm! was selected. This paper 
: presents a statement of the algorithm, boundary condition 
specifications, and the results of computational 
experiments compared with experimental data, j 

* H. Problem Statement j 

The three-dimensional PNS equation system for the 

* steady, subsonic turbulent flow of an isoenergetic fluid, to 

the principal scale of ordering 1 , is • 
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| With the basic governing, equation system thus 

idefined, construction of a suitable algorithm requires 
'addressing the ordering analysis and the fariieid 
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f The variables appearing in eqi|iations 1-3 have their 

• usual interpretation in,fluid mechanics, where superscript 
bar denotes conventional time averaging^, and the cross 

■ coupling of fluctuations in density and velocity have been 
! assumed i negligible. The tensor index summation 

convention i is implied, with xj aligned with the principal 

■ flow direction* and li < j < 3 and 2< (k,t)< 3- The 

'turbulence kinetic energy k is the trace of the Reynold* 
jitress tensor, C is the isotropic dissipation function, and Re 
lil the character l*lc Reynolds num&er. _-—--—' 
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For the reported results; the Reynolds stress tensor 
|field constitutive model of Biker et.ali 1 is employed, 
jlhe PNS ordering analysisindicates the extremum 
■significance of components of ufuj'is one order smaller 
-than unity. Simplifying the constitutive equation to this 
order yields 


I 


111. Finite Element Penalty Algorithm: 
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■The coefficients 

c Q . i 5 ° 

< k; in 


( 6 ) 


0.067, 0.56, 0.061} were used !ior the present study. 

. In the primitive form, equations 1-6 do not 
■ represent a well-posed initial-boundary value problem 
discretion lor the subsonic flow, multiple free-jet 
;problem class. As a consequence of the ordering, the 
scalar continuity equation governs both components of 
velocity in the transverse plane perpendicular to xj. As 
cited in the Introduction, various algorithm constructions 
•have been formulated to address this issue. The 
approach taken in the present analysis is to employ a 
finite element penalty function formulation, and: to 
append the order ( 6 } transverse plane momentum 
equations, 

; t 4 ( 5 k y * 517 [?*'"[] 4 1=7 • e !^ ' 0 Wi 


.to equation 3. Further, equation 3 is rearranged to the 
pressure Poisson equation. 


j Theoretical Statement 

* For the dependent variable set qUj)V {u,, k, c. p, 
jtljuj 3, equations 2, 3 * 7, k; 5, S; 6 ♦ 9 represent a weil- 
: posed, initial-boundary value problem description on th* 

. three dimensional domain >1 s R 2 i xf < {« , :i,C R 2 

: and x[ c C x j,x j}}. There is the additional fundamentali 
requirement that equation I be rigorously satisfied, since 

* It governs first order phenomena. The problem statement 
'i% completed with a polytropic equation o! state for 

determination of p for the isoenergetic flbw. j 

, As the consequence of the PNS reformulation, each, 
.of the first; six members of the set q are eligible for 
; constraint, on the boundary 3 R of R 2 , by a linear 
combination, of Dirichlet and Neumann boundary 
conditions. The first five of these members are alio 
.required specified as an initial-condition en the plane 
R 2 x x^. No boundary conditions are appropriate fbr the 

* algebraic equations governing 

The finite element penalty algorithm, for 
determination of the semi-discrete approximation q^(x^) 
to q(xj), is based on classical concepts for ciffereruial 
:constraints in the statement of variational ibsundary value 
■ problems. These concepts are extended to me very non¬ 
linear PNS problem class using a Calerkin, weighted* 
.residuals formulation. Deferring details 7 , t=e transverse 
plane domain R 2 is discretized into the union of non- 
.overlapping subdomains R 2 , wherein the ftmctionaJ form 
for the dependence of the semi-discrete •pproxitnaiton 
Tj h is assumed a priori specifiable. A convergent form is 
the cardinal basis { N k (x.)}, the members o! which are 
'polynomials on x^ complete to degree k_ Hence, the 
semirdiscrete approximation becomes tfbs union, ol 
elemental approximations, 


« ( *j> 1 «■»«/'» y; 

«* k ,U t ) 3 T {QK*.B k 


UO) 


01) 


where subscript and/or subscript e denotes perxaining to 
• the (finite element) domain Jl 2 . Ftsmer, {Ql} e 
’represents the values taken by q* at the: nodes of the 
domain, R 2 , and I <_ 1 < 12 is a tensor index denoting the 
.appropriate (nodal) vector dependent variable set of q^j 


Ifgv S }_! (u ) • l _ IjfL ♦ l— fo - 0 (S) ii With definition of'qh, equations 10-IT pesrmit direct 

^ k 5*^ [ k ij • evaluation ,of the semi-discrete approximation error LJC^hj 

^ ii ‘in each of the PNS governing differential equations. The 

* jibasic concept in the caltulus of a discreuxec variational 
Equation ft represents a quasi-linear elliptic ’ ‘ ** 

.boundary value problem, possessing complementary and! 

[particular solutions!. For the free-jei problem class; the 
; Complementary SDldtion is a homogeneous constant equal 
I to the farfield reference pressures The particular 
| solution to equation ft is thus obtained using 

'homogeneous Dirichlet boundary conditions. Upon 

|addition of equations 3 and 7* L(u k ) ♦ l_ 6 (u k ) represents a 
well-posed initial-boundary value problem for Gk, upon 
addition of the order (6 2 ( terms to the appropriate scalar 
^components of the Reynolds stress tensor, which are j 


| DAftilp LVl'KwCpi in |l»c LaiWUiUt 

'boundary value problfem is to render this error extremum 
!:in some norm. In the Calerkin we=gned-residua!s 
{extension of this concept, this is accomplished by 
requiring this error to be orthogonal tz> the space of 
'functions (iN k ) selected to define the sec*-discrete 
* approximation, i:e. 
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I The middle expression in equation 12 emphasizes that 
Jthe actual caJ cuius operations are performed on the 
■elemental domains Rjf. The resultant element (column) 
jmetrices are projected to the matrix structure ol the 
'global domain using the assembly operator 5 t , which is 
Jaimply matrix addition by rows 7 . t 

I Eqyatioo 12 defines the numerical solution 
JaJgorithm for the complete set with the exception el 
,‘the combined equations 3 and 7 for Mg., Here, the 
^definition of the extremum i must bi augmented 
^penalized) such that the continuity equation, is also 
satisfied! The functional form for this statement, which 
;i» an extension of the classical concept 7 , is 

EM .J 

» s/W V . UJl 

i *4^^™°* ; 

I 

The actual calculus operations defined in equation 13 are 
again , performed on an elemental basis and assembled, 
and X is an (arbitrary) parameter penalizing the 
statement of semi-discrete error orthogprsalixation for 

«V 

t: Equations 12-13 define the finite element penalty 

algorithm lor the PNS equation system. The theoretical 1 
arbitrariness remaining is solely the degree k of the 
cardinal'basts (Nj*),, spanning either three-sided or four- 
sided-element domains RZ, and the penalty parameter X 
and functional form for L(E*). However, equations 12413 
are definitions ol non-linear matrices, and resolution of 
the resultant problem definitions in linear algebra 
remains. For the semi-discrete approximations (u^.W h , 
C h ), equations 12-13 are matrix statements expressing 
the xj-ordinary derivative of the appropriate elements 
Ol (Ql) e , ^y*^‘o n U. A Taylor senes defines the 
matrix algebra statement lor the assembly of these 

elements of ( Ql), 1< lc 3, as. 


•where p is the iteration inddx at step xj,i,,and 


<«»$ ; ««»£i ♦ <«'»£}, 


M s n 


Equation 17 defines the fully-discrete approximation to 
the ^pendent variable set at the nodes of UR 2 * hence 
also q n (xj) throughout R 2 , see equation 10. Equation IS 
defines the (Newton) Jacobian of the non-linear algeoraic 
equation system, equation 13. 


Some Basic Decisions . I 

I Equation 13 delineates the basic decisions to be 
made regarding implementation of the penalty algorithm, 
into a computer code. In addition, for G > 0, ,a decision is 
required regarding approximate construction of the 
Newton algorithm Jacobian [ J],, equation 13, since it's 
size for the twelve dependent-van ad Id PNS statement is 
unwieldy on present computers. The decision on k, 
equation, 111 of course impacts considerably on the *»ize" 
of [31 

The multiple free jet analyses reported herein were 
Conducted using the GMC:3DPNS computer program .*’ 10 
Each of the basic decisions have been made ana evaluated 
lor this codt. The discretization of R 2 is cefined as the 
union of triangular cross-section finite eiem-ms spanned 
by the linear (k«l ) natural coordinate cardinal basis. The 
trapezoidal rule is empldypd for the integration 
algorithm, 0 * h in equation 14. The penafry, parameter 
X, following extensive numerical experimentation, is 
defined as the diagonal matrix. 


.{'F3> = {W iH - (QI> 4 • AxlQl)^ - 


M = C Ax R/ljJ +1 


In equation 1*, (QU j*| represents this ordinary 
derivative evaluated at some location on the interval 
Xj*| - Xj * Ax as defined by the parameter 0 < 6f 1. 

1 Equation 12, evaluated for, ph and u]ujV yields 
directly the appropriate column matrix statement { FI}.'« 
f{02. S < l< 12. Combined with equation 1$, the resultant 
linear "algebra statement of the finite element penalty 
algorithm for the PNS equation system, becomes 

, . |n[k. i; e. to. - wj <1}> 

Equation 13 Is a highly non-linear algebraic 
equation system, the character of which is largely, 
determined by the choice of the arbitrary solution 
parameter* k, X t 0, and'Ax. Equation 1,3 ddes not readily 
admit a useful Linearization, even for 0 ■ 0 which 
corresponds to explicit integration. Hence, the 
appropriate solution statement is the (Newton) matrix 
solution form. 


1 jfim]j P 
l ‘ ! 3s 


{Mil??! - - {rn p 


where C is a constant of order unity, and the elements (on 
the diagonal) of (“Ultf.j are (ODjLx* the nodal 
distribution of (Ql) computed at each iteration p at Xj,|. 

The functional form of the penalty term, equation 
;13, involves definition of,an auxiliary dependent■ variable 

♦ h . as 


, »:* f-h-hl . 

l(p 1 * X5J [P “jj 


The boundary conditions for 4 are a linear combination 
of homogeneous Dinchlet and Nfumam constraints, 
defined according to required flow porx?wty on various 
segments of ?R; Hence, Tjh is augmented for one 
additional entry, the solution ol 


, l 
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I The Galerkin weighted-residuals algorithm 
■statement for equation 21 is the algebraic equation 
•system <F13): « {0) t which is added to equation 13, 
.hence equations 1&-1S. Therefore, the explicit form of 
[the complete penalty term in equation 13 is. 




»{•<,> -h. 

T-i- Uo lilt 




• Regarding the Newton algorithm Jacobian, 
lequationi IS is replaced with two sparse matrices, 
.yielding a corresponding compromise on overall 
convergence rate while significantly reducing the matrix 
Yank; The initial-valued dependent variables {uty,k, c) 
? ar« sequentially solved as multiple ^ right side 
substitutions to equation! 16, using the Jacobian, 

[31 i] j*l* where 


mn - ?fn > 
IJ11] a SIyI7. 


The Poisson field variables {pf\ ♦ are solved 
sequentially as multiple right tide substitution s usi ng the 
pfr Jacobian [366]. The scalar components of u t uj are 
determined using an elemental averaging and assembly 
procedure, equivalent to solving equation 16 using [377] 
lor multiple right side substitutions. The algorithm 
timing utilises the sequence [311],, p66J , [377 ],,with 
updite ol the Jacobians occurring at every iteration. 
Detail! of the formation of these 3acobians is given m 
reference 7. 


IV. Results and Discussion 
Benchmark Tests 

The PNS finite element penalty algorithm, as 
operational , in the CMC:3DPNS computer code, is well 
documented for subsonic aerodynamic problem 
definitions involving semi*bourvded and fully bounded 
solution domains fl M*. Additional documentation is 
reported* 2 for *jn unbounded aerodynamic wake-type 
flow. Each of these includes detailed comparisons 
between algorithm prediction and experimental dita, 
including: the complete Reynolds stress tensor for two- 
and three-dimensional geometries. These comparisons 
'serve to quantitatively verify the appropriateness ol the 
■standard’: turbulence model constants ,Q* f equation 6, 
^which are used for the free-jet analyses. 

j As the typical case in three-dimensional, turbulent 
aerodynamic flbw prediction, quality, experimental data 
-to initiate solutions is usually non-existent. Hence, tests 
(were conducted to evaluate self-generation of initial 
Conditions for the free-jet conliguration. The banc 
'assumption was that the sole available initial data is the 
'dominant velocity scaJax component u|(«|,s.) on the 
PNS solution initiation plane. Figure 2 illustrates a two- 
dimensional sl^t jet test problem thus characterized by 
- u t* within and ex tenor to the jet. Since a 
non-zero background level of; dissipation function, isa 
computational I requirement (c-I appears throughout 
•’equations, 4-6), background levels ol both k° > 0 and 
*®>0 are defined to yield a background “turbulent 
Lviscosit)T_ levelk?/c _pi the . or der_ol_Jhe — 


i laminar viscosity v . Numerical experimentation 
.indicates:kO a OUO' 4 ) and C ° * W10* 9 ) yield aoequate 
. algorithm Stability,. At the knee of the prolife m ufl. 
Figure 2, the initial'Itvei ol turbulent viscosity is assumed 
; known, e.g.. 10 < v*/u < ID?- Since the nominal 
-extremum order of:Ic is 0<lO -2 ), this can be achieve* by 
Kiting k° • CKlU‘2)and t°= 0d0~ 7 ). 

Figure 3 summarizes the: PNS penalty algorithm 
prediction for the symmetric-half slot jet problem, for 
Vlp»er * 0.02 u? et i ^upper v? * WIO 2 !, on the 

span 0 < xj/Hjir 1L0, where H[ is the slot jet half-width. 
The initial condition for uj was interpolated as a step 
function on the nodes of the discretization of 12 , yielding, 
the spiked initial conditions lor k° and c® shewn as solidi 
lines. Tihe maximum levels of k and C increase by a 1 
factor of 2-4, by *i/Mf * 0.3, and thereafter decrease 
monotonically as the ennaheed level ol v t diffuses into 
the velocity defect region; The jel potential core is 
eroded by about 23*> by x|/Hji* 1.0. 

Figure 3b) summarizes the entrainment action 
predicted by the penalty algorithm. By definition, uj. 5 
.. 0, as shown by the solid line; Since the order delta uj 
momentum equation is homogeneous, the sole source for 
uj V O is the action ol the penalty constraint in equation 
13, The PNS solution indicates a large entrainment at 
I|/H{, = 0.3,, where u? = -0.1 in the farfield, which i 

progressively moderates as the solution proceeds further 
downfield. The boundary conditions lor; e for this 
geometry are 4 * 0 at the larheid and 3$/2*2- 0 on the 
symmetry line. 

The second basic test is the three-dimensional 
problem of a single jqt o! initially, Circular cross-section. 
For this case, Uj el s uf (k 33 m/s) and C?*tcT, or = O.i uft 
Figure 4 summarizes the PNS prediction ol transverse 
hall-plane velocity field at x|/Hj * 1.0, where Hj is the 
radius ol the initial jet, for laminar flow and for turbulent 
flow with v® * liO . Both solutions, exhibit the required 
symmetries and predict essentially radial entrainment. 
The plot length of each individual velocity vector is the 
measure ol the relative magnitude ol , scaled to the 
local Ipredicted extremum magnitude uy\i For the laminar 
flow prediction, u.J\ * 0:0317,, while 0^ * 0.044 for the 
turbulent case. thus, the measure ol magnitude of 
entrainment lor the turbulent case is approximately, 20 
times that for the laminar case, in qualitative agreement 
with expected behavior., The boundary conaitions lOf + 
for this case are f ■ 0 everywhere in the fariieid and 
^♦/2 x 2 ■ 0 on the symmetry plane; As in the two- 
dimensional test, = 0, andl the penalty term in; 

equation 13 is the sole causal mechanism initiating and 
maintaining the computed transverse plane velocity field. 


Multiple Free Jet PNS Prediction 

! The multiple jet case of main interest correspond* 
to a symmetric lour-jet geometry with the jets located in 
l,Close proximity and of small initial diameter. This 1 
configuration is verified, using experimental smoke flow 
visualization techniques, to rapidly induce a substantially 
.larg? transverse plane velocity field which efficiently 
pumps fluid, initially interior to the circvmierence ol the 
jets, into the exterior region. Figure 5*3 illustrates,the 
persistant unidirectional flow o( smoke obtained .with the 
multi-jet system inoperable. Figure 5b) vhows the rapid 
smoke dispersal 1 promoted by the multiple jet system 
operating at design conditions. 
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I Figure 4 Is * layout sketch of the multiple jet 
•device: indicating characteristic dimensions. Operating 
[on design, free air is induced to How down each vent 
'channel,, o! length 0.02 m, it a nominal velocity 
|Gj « 12„m/». The locator radius of the jets is R * 0.004m, 
land the initial hydraulic diameter of each jet is d; * 
*0.001 m. Ih the region interior to R, ui - 0j^2 u£;whiie 
in the essentially unbounded exterior region uj * u. The 
[characteristic Reynolds number for the jet Ho* ** 
,Re s I0 6 /m, and the vent channel walls are quite 
"‘rough. 

j The symmetry of this multiple jet geometry 
’permits the majority of PNS calculations to be 
performed on the symmetric quarter domain with 
boundary QABC, Figure 6 j Figure 7 graphs this domain 
'showing the nodal coordinates oi the basic M r 19 a 19 
‘non-uniform computational mesh. The lateral extent of 
the domain spans twice the locator radius R. Segments 
OA and 03 are symmetry planes, upon which the normal 
component of velocity vanishes and all other variables 
possess vanishing normal derivatives; Boundary segment 
ACB Is assumed sufficiently remote. Such that ail 
dependent variables have vanishing normal derivatives 
except ♦, which vanishes identically since the boundary 
U porous. 

For the bafic assessment, the sole specified initial 
conditions are ), and the levels of in the air 

jet and background flows. Since the PNS predicted 
secondary vortex flow field develops very rapidly; a 
relatively Ikrge (20%) background uj velocity field was 
used for efficiency. (The PNS algorithm stability was 
marginal using a 10% background llow specification. 
Small, but not significant differences were predicted in 
U a (x |, i j ) distributions at various stations downstream 
of the injector lice. The assumed initial level of v t 
exerts a much more significant influence.) Therefore, 
within the air jet, uf * 1.2 U w , where U. * 12. m/s, the 
design velocity. Everywhere exterior to the locator 
radius, uf I 0.2 V», and interior = (0.2 * 0.02) 
U a . Hence, the u^ velocity strain rate distribution on 
the initial surface plane was on-design. Further, witlun 
the air jet k° = 0.005, with t° defined such that 
v*/v r 35; Everywhere exterior to the jet, k° * OjOOOI 
and v°/v s 3. 

. Figure I is a composite of the multiple vortex 
secondary velocity field predicted by the PNS algorithm 
•t i|/R « 1,5, ie., 6 mm downstream from the plane of 
solution initiation. As in Figure 4, each velocity vector 
length is scaled to the extremum predicted level 1 
(u t /UU )max = uT * 0.067, The original location of 
[each air jel coincicfes with the clustering of large radial 
jvelocitiei, and the net action of the four-jet device is to 
induce a systemi of eight symmetrically disposed, 
’counter-rotating transverse vorrtex pairs. At this X) 
station, the maximum axial velocity component within 
|the initial jet region was uj « 0.33, hence the lc*_*i ratio 
[uTVo^ « 0j19 is quite substantial in comparison to the 
jPNS ordering analysis. 

| Quarter-plane plots of the PNS predicted evolution 
[of the transverse plane vortex flqwfieid in the near 
jvicinity of the locator radius and air jft are shown in 
’Figure 9. A substantial field exists by xj/R * 0.3; which 
: U 2 mm downstream from the initial condition plane, 
.although the vortex pattern has not matured. However,, 
.by xi|/R ■ 1.0, Figure 9b), the counterrotating pattern is 
^clearly evident and « 0.051 ha* remained essentially 
-constant. By xj/R . 1.3, Figure 9C, the visual 

-appearance of the u t distribution U nominally 
Ijjnchangpd, although the overall strength has increaisedl 
J>y.aboja.l3>a-Wnce,uY «_0.D4£._The visual appearance_ 


of the double counter-rotating vortex pattern shown in 
Figure 9c) persists essentially unchanged if or. distances up 
to xi/R « HI ddwnatre*m of the initial condition planer 

A species continuity equation, similar in appearance 
.to equation 7, was added to the PNS.system, to permit . 
tracking of the distribution ol flbid initially interior to 
the locator radius R. Figure 10a shows the Y P « 103% 
Smoke distribution at the solution initiation pline; the 
four troughs correspond to the location of the four air 
jets. Diffusion processes dominate on 0 < xyJ R < l. 
Figure 10b), and the extremum: smoke concentratibn”Y ni 
remains 100%. By t\IR * 2.0, Y, m « 97*%,* at i,/R,r 3.3 
Y, m = 84%, and the first norwerD level Y e occurs at the 
boundary 3R of the PNS domain at xj/R * 5.0, Figure 
10c), where Y m; « 63%. By xj/R = 1II.0. Y m has decreased 
to 30% and Y® * 23%, Figyre lOd). Hence* tne counter¬ 
rotating vortex system has homogenized the initial smoke 
level within a distance of about 44mmj Of greatest 
significance, on }|R at xj/R = 11.0, 3Y/a:X| • Aj, > 0 

confirms that the material transport is dominated by 
convection, in qualitative: agreement:with the smoke flow 
visualization experimental data. 

Figure 11 summarizes the PNS computed 
distribution of turbulent kinetic energy k(xy ) on 
0 < xij/R < 1.5, The initial condition is kO = GiC05; 
Figure 1 la£ within the jets and:k<> * 0.0001 elsewhere. 
Due to the rapid decay of the air jets; the extremum level' 
of k it x|/R s 0.12 has increased sharply to k m = 0.023, 
Figure U:b). This has further, increased to = O.C35 it 
xji/R - 0.25, Figure Mc^and the: crater s-hape: of each 
profile is clearly evident. Diffusional processes 

eventually smooth the distributions, Figure lid), and the 
extremum levels has been reached,, km * 0.039. 
Continuing downstream, the distributions of k become 
further homogenized at nominally the same maximum 
level. 

A series of computational experiments were 
conducted to ascertain importance ol assumed initial 
turbulence level on the PNS prediction. For one test, ko 
was halved to k° = 0;0025 and c° held constant, yielding 
VtAi » 9. A second test was conducted with the air jets 
assumed laminar. Figure 12 summarizes these predictions 
in terms of decay of the jet velocity extremum, ,i^(xj), 
•nd the extremum predicted vortex velocity component 
For the laminar flow case, the air jet extremum 
does not decrease at all on 0 < xi/R<_ 1.5, and the induced 
.transverse plane component “hover® about u7: 0.0 L. In 
distinction, halving the initial turbulence level simply 
'displaces the jet decay curve to the right, a distance oi 
!Akii/R ti 0.25, and yieldi an extremum transverse 
;component that is only modestly smaller by about Au'f : 

10.00). In comparison to the experimental dita. these 
’predictions confirm the importance of the vent channel 
roughness in promoting the desired action of the multi-jet 

• system. In the same sense, the assumed initial level for 
;k° dies not appear critical, in terms ol the PNS soHihon, 
'predicting a qualitatively valid solution, provided vt '* 
[sufficiently large to permit self-generation of solutions to 
[the k-c equation system. 

* Additional experimental data confirms that 
sufficient protrusion ol the locator radius surface beypnd 
the multi-jet injector lace can markedly reAjcc the level 
of visually apparent Swirling motion. This geometry 
modification was modeled by assuming (he locator radius 
R, Figure 7, was a sol.d surface for a speciHed distance 
Ax(/R downstream. Hence, the initial PNS solution 
domain is contained within the radius surface, and *he 
boundary conditions on this surface are q-* ■ ( 0!)* After 

..marching the prescribed distance, the PNV solution is 
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| stopped, The locator radius surface is computationally 
I removed and the remainder of the solution domain OABC 
■ added i and the PNS solution restarted using the current 

I solution as initial conditions. j 

Figure 13 summarizes the resulting PNS prediction 
' of u Since the solid locator radius surface is a 

!*ink Joe momentum, one expects that a proportionally 
lirger level 1 of transverse plane velocity will result. 

J Figure 13a) confirms this for i 1 mm extension. The 
nominally radial flow from the jet region is evident, it 
does not penetrate the radius boundary,,and uT = 0.112 
is more than .double the extremum predicted without the 
extension. Figure 13b) shows the u t distribution 
computed at xj/R ■ 0 j 75 with the solid surface 
appended! A weak vortex pair is barely discernable, and 
u"V* 0.032 has decayed by a factor of three. Figure 13c) 
Shows the u t distribution computed at A*j/R * 1.3, 
downstream of the limm extension, which used the data 
of Figure 13a) as initial conditions. Comparing to the 
on-design solution. Figure 9c), a single vortex pair of 
much weaker strength has replaced the basic prediction. 
Even though the extremum of u™ ■ 0.055 is only 
modestly decreased, the vortex pattern of Figure 13c) 
would be much less effective in promoting convection of 
. material from within the interior region. 

i 

For a final I computational experiment, the 
available smoke flow visualization data indicate that the 
multiple jet device efficiency is only modestly aitered by 
partial blocking of one of the initial jets. Since this 
constraint destroys the inherent geometric symmetries, 
the PNS solution domain must now encompass a region 
that is a factor of four larger. Figure 14 graphs the 
transverse velocity distribution Dg at xj/R = l.J on an 
M * 19 x 19 mesh that is twice as Coarse as the base 
discretization; Figure 7. However, comparing Figure S 
and lk confirms that the coarse grid full solution has 
captured the essential! eight-vortex pair structure, and 
GT « 0.05a U. within 20% of the liner grid extremum, 
Figur e 9 c). 

Figure 13 graphs the PNS predicted transverse 
plane velocity distribution ui. at xi/,R c 1L3 obtained 
with the lower right jet completely shut off; A through- 
flow has resulted over the occluded jet, and the 
extremum level'of * 0.12 is about 50% larger than 
the: design. configuration. Figure 16 compares the 
computed distributions of smoke density at x^/R > 6.0 
_ for both coarse grid solutions: Figure 16a) compares 
visually with the finer grid solution. Figure 10c), and Y m , 
:» 53% lies on the interpolation of trajectory extremum, 
.The loss action of the occluded vent is deiriy evident in 
.Figure 16b), although Y m * 39% is.only 10% larger than 
* the on - design , solution. These PNS predictions thus 
agree qualitatively with the field data, and further 
permit a quantitative comparison measure of the action 
•■of design modif ications. 


I V. Condusiont 

! The finite element penalty numerical solution 
j algorithm for the three-dimensional parabolic Naviee- 
iStokes,equations, for subsonic turbulent flows, has been 
♦applied to prediction: of secondary vortex flowfields 
induced by, multiple free-jets issuing in close proximity. 
,The combined action of rapid jet decay and discreteness 
of the multiple jet iconfiguration has been quantitatively 
t assessed I regarding resultant entrainment and secondary 
‘vortex structures. No detailed experimental! 

measurements, are available for comparison of the 
: numerical predictions. However, interpretation of 
i_ir>eipensj*ely_acquired_vtdeo-graphi c s moke flbw • 


•visualization data has provided a basis for qualitative 
; comparisons. Further, a range of comprjrationali 
» experiments were conducted; involving parameter 

; variations, that permitted further qualitative 

comparisons. The availability of the computer program 
i embodiment of 1 the theory, as a computational laboratory, 

1 has been verified of utility and of sufficient versatility to 
permit a range of experiments. It is fair to assume that 
. this constitutes one small step toward* the eventual 
realization of computational fluid dynamics as a 
} diagnostic practice. 
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f igure 3. Penalty Algorithm Solution of Symmetric- 
Hail Subsonic Slot JlJ, uj = 30 m/s, O 1.0, a) 

Axial Mean Velocity u\ t U) Transverse .VTean Velocity »2» 
3) Turbulent Kinetic Energy k, d) Isotropic Dissipation 
Function g. 
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figure 4. Penalty Algorithm Solution For Transverse 
Mean Velocity up Distribution, Symmetric-Half Circular 
free Jet, ui = iS m/s, xj/Dh = 1.0, a) Laminar Jet, u»J - 
0.0017, b) Turbulent Jet, uj = 0.044. 
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Figure 5. Smoke Flow Visualization ol Multiplc-3et 
Configuration, a) jet Blows Oil, b) Jet Flows On. 




Figure 6. Engineering Layout of Multiple 3et 
Geometry, a) Plan-View, b) End View, Dimensions in mm. 
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Figure 7. Penalty Algorithm Solution Domain For 
Symmetric Quarter Plane Prediction of Four Multifile 
let Configuration Including Nodal Coordinate 
Distribution of M = 19 x 19 Discretization. 
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Figure 8. Composite Penalty Algorithm Prediction of Transverse Plane Velocity T3^ Distribution, Four 
Geometry, ui s 12 m/s, xj/R = 1.5. 
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Nearfield Transverse Plane Velocity 
ih Distributions, Penalty Algorithm Prediction, uj = 12 
m/s, a) xj/R = 0.5, u^ 1 - 0.053, b) xi/R = 1.0, u'f 1 - 0.058, 
c) xj/R = 1.5, u r J> = 0.067. 
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